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Viscous stabilization of 2D drainage displacements with trapping 
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We investigate the stabilization mechanisms due to viscous forces in the invasion front during 
drainage displacement in two-dimensional porous media using a network simulator. We find that 
in horizontal displacement the capillary pressure difference between two different points along the 
front varies almost linearly as function of height separation in the direction of the displacement. 
The numerical result supports arguments taking into account the loopless displacement pattern 
where nonwetting fluid flow in separate strands (paths). As a consequence, we show that existing 
theories developed for viscous stabilization, are not compatible with drainage when loopless strands 
dominate the displacement process. 
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Immiscible displacement of one fluid by another fluid 
in porous media generates front structures and pat- 
terns ranging from compact to ramified and fractal [1-3]. 
When a nonwetting fluid displaces a wetting fluid 
(drainage) at low injection rate, the nonwetting fluid gen- 
erates a pattern of fractal dimension equal to the cluster 
formed by invasion percolation [4]. The displacement is 
controlled solely by the capillary pressure, that is the 
pressure difference between the two fluids across a pore 
meniscus. At high injection rate and when the viscosity 
of the nonwetting fluid is higher or equal to the viscosity 
of the wetting fluid, the width of the displacement front 
stabilizes and a more compact pattern is generated [2,5] 

The purpose of the present letter is to investigate 
the stabilization mechanisms of the front due to viscous 
forces. To study the stabilization mechanisms we consider 
two-dimensional (2D) horizontal drainage at different in- 
jection rates. Since the displacement is performed within 
the plane we neglect gravity. We present simulations 
where we have calculated the capillary pressure difference 
AP C between two different pore menisci along the front 
separated a height Ah in the direction of the displace- 
ment [Fig. 1(a)]. The simulations are based on a network 
model that properly describes the dynamics of the fluid- 
fluid displacement as well as the capillary and viscous 
pressure buildup [6,7]. Simulations show that for a wide 
range of injection rates and different fluid viscosities AP C 
varies almost linearly with Ah (Figs. 2 and 3). Assuming 
a power law behavior AP C oc Ah K we find k = 1.0 ± 0.1. 
This is a surprising result because the viscous force field 
that stabilizes the front, is non homogeneous due to trap- 
ping of wetting fluid behind the front and to the fractal 
behavior of the front structure. 

Based on the observation that the displacement struc- 
tures are characterized by loopless strands of nonwetting 
fluid [Fig. 1(a)], we also present arguments being sup- 
ported by our numerical findings. We conjecture that the 
arguments might affect the behavior of the front width 
w s as function of the capillary number C a . Here C a de- 




FIG. 1. (a) The displacement pattern of a simulation on 
a lattice of 25 x 35 nodes. Nonwetting fluid (dark grey and 
black) is injected from below and wetting fluid (light grey) 
flows out along the top row. In the figure, AP C is the capillary 
pressure difference between a meniscus at the bottom filled 
dot and a meniscus at the topmost filled dot, separated a 
height Ah in the direction of the displacement. The black 
tubes indicate strands containing no loops, where nonwetting 
fluid flow. The dark grey tubes connecting to the strands, 
are dead ends where nonwetting fluid cannot flow because of 
trapped wetting fluid, (b) A tube filled with wetting fluid and 
surrounded on both sides by nonwetting fluid is trapped. 



notes the ratio between viscous and capillary forces and 
in the following C a = Q[i nw /Y,j, where Q is the injection 
rate, S is the cross section of the inlet, and \x nw is the 
viscosity of the nonwetting phase. 

In the literature [8-11] there has been suggested 
slightly different scaling behavior of w s as function of 
C a and a general consensus has not yet been reached. 
However, none of them consider the evidence observed 
here that the displacement patterns are loopless and that 
nonwetting fluid only flows in strands to displace wetting 
fluid. As a consequence, we show that earlier proposed 
theories [8-11] can not be used to describe drainage when 



loopless nonwetting strands dominate the displacements. 

Before we present the numerical results and the theo- 
retical evidence, we briefly introduce the network model. 
The model porous medium consists of a square lattice of 
cylindrical tubes oriented at 45° to the longest side of the 
lattice [Fig. 1(a)]. Four tubes meet at each intersection 
where we put a node having no volume. The disorder is 
introduced by (1) assigning the tubes a radius r chosen 
at random inside a defined interval or (2) moving the in- 
tersections a randomly chosen distance away from their 
initial positions. In (1) all tubes have equal length d but 
different r. (2) results in a distorted square lattice giving 
the tubes different lengths. Here r = d/2a where a is 
the aspect ratio between the tube length and its radius. 

The tubes are initially filled with a wetting fluid of vis- 
cosity fj, w and a nonwetting fluid of viscosity (i nw > jjl w , 
is injected at constant injection rate Q along the bot- 
tom row (inlet). The viscosity ratio M is defined as 
M = Unw/Hw The wetting fluid is displaced and 
flows out along the top row (outlet). There are peri- 
odic boundary conditions in the orthogonal direction. 
The fluids are assumed immiscible, hence an interface 
(a meniscus) is located where the fluids meet in the 
tubes. The capillary pressure p c of a meniscus is given 
by Pc — (27/V) [1 — cos(2ttx / d)] . The first term is Young- 
Laplace law for a cylindrical tube when perfect wetting is 
assumed and in the second term x is the position of the 
meniscus in the tube (0 < x < d). Thus, with respect to 
the capillary pressure we treat the tubes as if they were 
hourglass shaped with effective radii following a smooth 
function. By letting p c vary as above, we include the ef- 
fect of local readjustments of the menisci at pore level [6] 
which is important for the description of burst dynam- 
ics [12]. The detailed modeling of p c costs computation 
time, but is necessary in order to properly simulate the 
capillary pressure behavior along the front. 

The volume flux qij through a tube between the zth 
and the jth node is given by Washburn equation [13]: 
Qij = -((Tijkij I 'Hij){Pj - Pi —p c ,ij)/dij. Here % is the 
permeability of the tube, Uij is the average cross section 
of the tube, pi and pj is the pressures at node i and j re- 
spectively, and p c ,ij is the sum of the capillary pressures 
of the menisci inside the tube. A tube partially filled with 
both liquids, is allowed to contain one or two menisci. 
Furthermore, /i,j denotes the effective viscosity given by 
the sum of the volume fractions of each fluid inside the 
tube multiplied by their respective viscosities. Inserting 
the above equation for q^ into Kirchhoff equations at 
every node (volume flux conservation), Y] . q^ — 0, con- 
stitutes a set of linear equations which are to be solved 
for pi. The set of equations is solved by using the Conju- 
gate Gradient method with the constraint that Q is held 
fixed. See Refs. [6,7] for details on the numerical scheme 
updating the menisci and solving pi . 

The front between the two phases is detected by run- 
ning a Hoshcn-Kopelman algorithm [14] on the lattice. 
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FIG. 2. AP C as function of Ah for high and intermediate 
C a with M = 100 on a lattice of 25 x 35 nodes, and for low 
C a with M — 1 on a lattice of 40 x 60 nodes. 



The front width is defined as the standard deviation of 
the distances between each meniscus along the front and 
the average front position in the direction of the displace- 
ment. AP C as function of Ah is calculated by taking 
the mean of the capillary pressure differences between 
all pairs of menisci separated a height Ah along the 
front. The capillary pressure difference between a pair 
of menisci is calculated by taking the capillary pressure 
of the meniscus closest to the inlet minus the capillary 
pressure of the meniscus closest to the outlet [Fig. 1(a)]. 
Figure 2 shows AP C as function of Ah for simulations 
performed at three different C a 's with M — 100 or 1. The 
simulations with M = 100 were performed on a 25 x 35 
nodes lattice with fj, nw = 10 P, /i w = 0.10 P, and 7 = 30 
dyn/cm. The disorder was introduced by choosing the 



tube radii at random in the interval 0.05<i < 



< d. 



The tube length was d = 0.1 cm. The simulations with 
M = 1 were performed on a distorted lattice of 40 x 60 
nodes where 0.02 cm < rfy < 0.18 cm and r^ = dij/2a 
with a = 1.25. Here fx nw = /i w = 0.5 P. To obtain 
reliable average quantities we did 10-30 simulations at 
each C a with different sets of random r^ or dij . 

From Fig. 2 we observe that AP C increases roughly lin- 
early as function of Ah. At lowest C a no clear stabiliza- 
tion of the front was observed due to the finite size of the 
system. At higher C a the viscous gradient stabilizes the 
front. The gradient causes the capillary pressure of the 
menisci closest to the inlet to exceed the capillary pres- 
sure of the menisci lying in the uppermost part. Thus, 
the menisci closest to the inlet will more easily penetrate 
a narrow tube compared to menisci further down stream. 
This will eventually stabilize the front. 

To save computation time and thereby be able to study 
AP C on larger lattices in the small C a regime, we have 
generated bond invasion percolation (IP) patterns with 
trapping on lattices of 200 x 300 nodes. The IP pat- 



terns were generated on the bonds in a square lattice 
with the bonds oriented diagonally at 45°. Hence, the 
bonds correspond to the tubes in our network model. 
Each bond was assigned a random number fij in the 
interval [0, 1]. A small stabilizing gradient g = 0.05 was 
applied, giving an occupation threshold tij of every bond: 
tij = fij + ghij [8,15]. Here hij denotes the height of the 
bond above the bottom row. The occupation of bonds 
started at the bottom row, and the next bond to be occu- 
pied was always the bond with the lowest threshold value 
from the set of empty bonds along the invasion front. 
The generated IP patterns are similar to the site-bond 
IP patterns in [16] and we assume they are statistical 
equal to structures that would have been obtained in a 
corresponding complete displacement simulation. 

When the IP patterns became well developed with 
trapped (wetting) clusters of sizes between the bond 
length and the front width, the tubes in our network 
model were filled with nonwetting and wetting fluid ac- 
cording to occupied and empty bonds in the IP lattice. 
Moreover, the radii r^ of the tubes were mapped to the 
random numbers fij of the bonds as r%j — [0.05 + 0.95(1 — 
fij)]d. Thus, 0.05d < r^ < d and we set the tube length 
d = 0.1 cm. Note that r^ is mapped to 1 — fij because in 
our IP algorithm the next bond to be invaded is the one 
with the lowest threshold value, opposite to the network 
model, where the widest tube will be invaded first. 

After the initiation of the tube network was completed, 
the network model was started and the simulations were 
run a limited number of time steps before it was stopped. 
The number of time steps where chosen sufficiently large 
to let the menisci along the front adjust according to the 
viscous pressure set up by the injection rate. 

Totally, we generated four IP patterns with different 
sets of fij and every pattern was loaded into the network 
model. The result of the calculated AP C versus Ah is 
shown in Fig. 3 for C a = 9.5 x 10~ 5 and M = 100. If we 
assume a power law AP C oc Ah K , we find k = 1.0 ± 0.1. 
The slope of the straight line in Fig. 3 is 1.0. We have 
also calculated AP C for C a = 2 x 10~ 6 with M = 1 and 
M = 100 by using one of the generated IP patterns. The 
result of those simulations is consistent with Fig. 3. 

Wilkinson [8] was the first to use percolation theory 
to deduce a power law between w s and C a when only 
viscous forces stabilize the front. In 3D, where trapping 
of wetting fluid is assumed to be of little importance, he 
suggested w s oc C a ~ a and a = v/(l+t — (3 + v). Here t is 
the conductivity exponent and /? is the order parameter 
exponent in percolation. Blunt et al. [10] used a similar 
approach, however, they found a = v j(\ + t + v) in 3D. 
This is identical to the result of Lenormand [9] discussing 
limits of fractal patterns between capillary fingering and 
stable displacement in 2D porous media. Blunt et al. also 
deduced a scaling relation for the pressure drop AP nw 
across a height difference Ah in the nonwetting phase of 
the front and found AP nw oc Ah t ' v+1 . Later on, Xu et al. 
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FIG. 3. log 10 (AP c ) as function of log 10 (A/i) for drainage 
simulations initiated on IP patterns at C a — 9.5 x 1CP 5 and 
M — 100. The result is averaged over four different runs and 
the error bars denote the standard error in the mean. The 
slope of the solid line is 1.0. 



[11] used the arguments of Gouyet et al. [17] and Wilkin- 
son [8] to show that AP nw oc ^h t,v+dE ~ 1 ~' 3/v j where d ^ 
is the Euclidean dimension of the space in which the front 
is embedded. They also argued that AP C = AP nw — AP W 
where AP W denoting the pressure drop in the wetting 
phase of the front, is linearly dependent on Aft, due to 
the compact displaced fluid [see Fig. 1(a)]. Thus, the re- 
sult of Xu et al. would in 2D predict AP C oc Ah 19 where 
we have used t — 1.3, v = 4/3, j3 = 5/36, and c^e = 2. 
Our simulations give AP C oc Ah K and k = 1.0 ± 0.1. Be- 
low we present an alternative view on the displacement 
pattern from those first suggested by Wilkinson. The 
alternative view is based upon the looplcss nonwetting 
strands and is supported by our numerical result. 

The simulated displacement patterns show that the 
nonwetting fluid contains no closed loops [Fig. 1(a)] be- 
cause wetting fluid may be trapped in single tubes, due 
to volume conservation [Fig. 1(b)]. Because of fluid trap- 
ping in single tubes, the invading fluid flows in separate 
strands that cannot coalesce. We note that the definition 
in in Fig. 1(b) can be easily generalized to 3D [18], since 
increasing the coordination number of the lattice does 
not change the trapping rule. Therefore, we expect loop- 
less patterns to develop in 3D lattices and our arguments 
that we present below should apply there too. We also 
note that trapping of wetting fluid is more difficult in real 
porous media due to a more complex topology of pores 
and throats there. Loopless IP patterns have earlier been 
observed in Rcfs. [16,19,20]. 

From Fig. 1 (a) we may separate the displacement pat- 
tern into two parts: one consisting of the frontal region 
continuously covering new tubes, and the other consist- 
ing of the more static structure behind the front. The 
frontal region is supplied by nonwetting fluid through 
strands connecting the frontal region to the inlet. When 



the strands approach the frontal region they are more 
likely to split. Since we are dealing with a square lat- 
tice, a splitting strand may create either two or three 
new strands. As the strands proceed further into the 
frontal region they split again and again and eventually 
they cover the frontal region completely [see Fig. 1(a)]. 

On IP patterns without loops [16,18,20] the length I 
of the minimum path between two points separated an 
Euclidean distance R scales like I oc R Ds where D s is 
the fractal dimension of the shortest path. We assume 
that the displacement pattern of the frontal region for 
length less than the correlation length (in our case w s ) is 
statistically equal to IP patterns in [16]. Therefore, the 
length of a strand in the frontal region is proportional to 
Ah Ds when Ah is less than w s . If we assume that on 
the average every tube in the lattice has same mobility 
(kijf/Jij), this causes the fluid pressure within a single 
strand to drop like Ah Ds as long as the strand does not 
split. When the strand splits volume conservation causes 
the volume fluxes through the new strands to be less than 
the flux in the strand before it splits. Hence, following a 
path where strands split will cause the pressure to drop 
as Ah K where n < D s . 

From the above arguments we conclude that the pres- 
sure drop AP nw , in the nonwetting phase of the frontal 
region (that is the strands) should scale as AP nw oc Ah K 
where k < D s . In 2D two different values for D s have 
been reported: D s = 1.22 [18,20] and D s = 1.14 [16]. 
Both values are consistent with our simulations finding 
K= 1.0 ±0.1. 

The evidence that k ~ 1.0 may influence the scaling of 
w s as function of C a . At low C a simulations show that 
AP C oc C a Ah 10 [21]. Here AP C denotes the capillary 
pressure difference when the front is stationary. That 
means, AP C excludes situations where nonwetting fluid 
rapidly invades new tubes due to local instabilities. At 
sufficiently low C' a the displacement can be mapped to 
percolation giving AP C oc / — f c oc S,^ 1 ^ [8,15,17]. Here 
/ is the occupation probability of the bonds, f c is the per- 
colation threshold, and £ oc w s is the correlation length. 
By combining the above relations, we obtain w s oc C a a 
where a = i//(l + vk). In 2D v = 4/3 and inserting 
k = 1.0 gives a « 0.57. At high C a we expect a crossover 
to another type of behavior since it is not clear if the 
mapping to percolation [8,15,17] is valid there. We note 
that Wilkinson's result [8] gives a « 0.38 in 2D. 

In summary we conclude that AP C oc Ah K where our 
simulations gives k = 1.0 ± 0.1. By describing the dis- 
placement structure in terms of loopless strands [16,20] 
we have argued that k < D s , where D s is the fractal di- 
mension of the shortest path between two points on IP 
patterns without loops. In 2D two values of D s has been 
reported (1.14 [16] and 1.22 [18,20]) and both arc con- 
sistent with our numerical result k ~ 1.0. We conclude 
that earlier suggested theories [8-11] are not compatible 
in situations where a loopless pattern with nonwetting 



strands dominate the displacement. We have also shown 
that a in w s oc C a ~", may be influenced by the evidence 
that k < D s . Work is in progress to investigate our ar- 
guments in 3D and the effect of loops on n. 
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